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1 Abstract 

2 Purpose: Regional lung volume change as a function of lung inflation serves as an 

3 index of parenchymal and airway status as well as an index of regional ventilation and 

4 can be used to detect pathologic changes over time. In this article, we propose a new 

5 regional measure of lung mechanics — the specific air vohime change by corrected 

6 Jacobian. We compare this new measure, along with two existing registration-based 

7 measures of lung ventilation, to a regional ventilation measurement derived from 

8 xenon-CT (Xe-CT) imaging. 

9 Methods: 4DCT and Xe-CT data sets from four adult sheep are used in this study. 

10 Nonlinear, 3D image registration is applied to register an image acquired near end 

11 inspiration to an image acquired near end expiration. Approximately 200 anno- 

12 tated anatomical points are used as landmarks to evaluate registration accuracy. 

13 Three different registration-based measures of regional lung mechanics are derived 

14 and compared: the specific air volume change calculated from the Jacobian (SAJ); 

15 the specific air volume change calculated by the corrected Jacobian (SAC J); and the 

16 specific air volume change by intensity change (SAI). We show that the commonly- 

17 used SAI measure can be derived from the direct SAJ measure by using the air-tissue 

18 mixture model and assuming there is no tissue volume change between the end in- 

19 spiration and end expiration data sets. All three ventilation measures are evaluated 

20 by comparing to Xe-CT estimates of regional ventilation. 

21 Results: After registration, the mean registration error is on the order of 1 mm. 

22 For cubical ROIs in cubes with size 20 mm x 20 mm x 20 mm, the SAJ and SACJ 

23 measures show significantly higher correlation (linear regression, average = 0.75 

24 and = 0.82) with the Xe-CT based measure of specific ventilation (sV) than the 

25 SAI measure. For ROIs in slabs along the ventral-dorsal vertical direction with size 

26 of 150 mm x 8 mm x 40 mm, the SAJ, SACJ, and SAI all show high correlation 

27 (linear regression, average = 0.88, = 0.92 and = 0.87) with the Xe-CT based 

28 sV without significant differences when comparing between the three methods. We 

29 demonstrate a linear relationship between the difference of specific air volume change 

30 (DSA) and difference of tissue volume (DT) in all four animals (linear regression, 

31 average = 0.86). 

32 Conclusion: Given a deformation field by an image registration algorithm, signifi- 

33 cant differences between the SAJ, SACJ, and SAI measures were found at a regional 

34 level compared to the Xe-CT sV in four sheep that were studied. The SACJ in- 

35 troduced here, provides better correlations with Xe-CT based sV than the SAJ and 

36 SAI measures, thus providing an improved surrogate for regional ventilation. 



37 Keywords: image registration, ventilation, lung function, tissue function, pulmonary 
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1 1 Introduction 



2 Regional ventilation is the term used to characterize the volume of fresh gas per unit time 

3 that enters or exits the lung at the acinar (gas exchange) level. Disruption of regional 

4 ventilation can reflect alterations to airways (physiological or pathological), alterations in 

5 parenchymal mechanics, changes to the muscles of respiration, body posture effects and 

6 inhaled gas properties. Thus, measures of regional lung mechanics can serve as a sensitive 

7 test of the status of the respiratory system and should be considerably more sensitive and 

8 informative than global pulmonary function test. Recent advances in multi-detector- row 

9 CT (MDCT), 4DCT respiratory gating methods, and image processing techniques enable 

10 us to study pulmonary function at the regional level with high resolution anatomical 

11 information compared to other methods. MDCT can be used to acquire multiple static 

12 breath-hold CT images of the lung taken at different lung volumes, or 4DCT images of the 

13 lung acquired during spiral scanning using a low pitch and retrospectively reconstructed at 

14 different respiratory phases with proper respiratory control [Hill [3]. Image registration can 

15 be applied to these data to estimate a deformation field that transforms the lung from one 

16 volume configuration to the other. This deformation field can be analyzed to estimate local 

17 lung tissue expansion, calculate voxel-by-voxel intensity change, and make biomechanical 

18 measurements. When combined with image segmentation algorithms [H O [6l [7] , functional 

19 and biomechanical measurements can be reported on a lung, lobe, and sublobar basis, and 

20 can be used to interpret regional lung function relative to specific segments of bronchial 

21 tree. Such measurements of pulmonary function have proven useful as a planning tool 

22 during RT planning [HI [9] and may be useful for tracking the progression of toxicity to 

23 nearby normal tissue during RT and can be used to evaluate the effectiveness of a treatment 

24 post-therapy [TO] . 

25 Early studies using CT to study regional air volume changes have proved to enhance 

26 our understanding of normal lung function. Several groups have proposed methods that 

27 couple image registration and CT imaging to study regional lung function. Guerrero et 

28 al. have used optical-flow registration to compute lung ventilation from 4DCT [111 112] 

29 with an intensity-based ventilation measure. Christensen et al. used image registration 

30 to match images across cine-CT sequences and estimate rates of local tissue expansion 

31 and contraction pLSj using a Jacobian-based ventilation measure. While they were able to 

32 show that their accumulated measurements matched well with the global measurements, 

33 they were not able to compare the registration-based measurements to local measures 

34 of regional tissue ventilation. Recently, Castillo et al. compared the intensity-based and 

35 Jacobian-based calculations of ventilation from 4DCT with the ventilation from ^^^Tc- 

36 labeled aerosol SPECT/CT [H]. A statistically higher correlation to the SPECT/CT 
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1 based ventilation was found for intensity-based based calculation over the Jacobian-based 

2 calculation. However, the comparison of the two techniques was based on the Dice sim- 

3 ilarity coefficient between the thresholded masks within 20% variation from the 4DCT 

4 and from SPECT/CT. Though their experiment is novel and important, since the average 

5 mask size is about 490.5 mL (with average subject exhale volume 2452.7 mL, and 5 sub- 

6 masks per subject), the comparison is more global than regional. In addition, as shown in 

7 Section 12. 4[ both the intensity-based and Jacobian-based ventilation measures are based 

8 on the assumption that regional lung volume change is due solely to air content change, 

9 which may not always be a valid assumption. Other factors, such as blood volume change, 

10 may also introduce the regional lung volume change. 

11 The physiologic significance of these registration-based measures of respiratory function 

12 can be established by comparing to more conventional measurements, such as nuclear 

13 medicine or contrast wash- in/ wash-out studies with CT or MR. Xenon-enhanced CT (Xe- 

14 CT) measures regional ventilation by observing the gas wash-in or wash-out rate on serial 

15 CT images [151 IIHl EZ] Xe-CT imaging has the advantage of high temporal resolution and 

16 spatial resolution and refiects a measure of fresh gas delivery to the gas exchange units 

17 of the lung. Although the dynamic Xe-CT method is limited in Z-axis coverage, requires 
IB expensive Xe gas, and is technically challenging, it serves as the gold standard of regional 

19 ventilation and can be used to compare with registration-based measures of regional lung 

20 function in animal studies for validation purposes. 

21 This paper describes three measures to estimate regional ventilation from image registra- 

22 tion of CT images: the specific air volume change calculated from the Jacobian (SAJ); the 

23 specific air volume change calculated by the corrected Jacobian (SAC J); and the specific 

24 air volume change by intensity change (SAl). We show that the SAI ventilation measure 

25 can be derived from the SAJ measure by making the assumption that there is no tissue 

26 volume change between registration volumes. We evaluate these three measures by com- 

27 paring them with a Xe-CT measure of ventilation in a regional basis (20 mm x 20 mm 
2B X 20 mm cube, or 8 mL). Among these three registration based measures, we show that 

29 the corrected Jacobian-based measure, SACJ, has the best correlation with the Xe-CT 

30 derived measure of specific ventilation. 
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1 2 Material and methods 



2 2.1 Method Overview 

3 Our goal is to validate and compare the measures used to estimate regional lung ventilation 

4 from image registration by comparing them to Xe-CT estimated ventilation. Figure [T] 

5 shows a block diagram of the entire process. Two types of data were acquired for each 

6 animal: a 4DCT scan and a Xe-CT scan. In order to make our comparisons under the 

7 same physiological conditions, each animal was scanned and mechanically ventilated with 

8 the same respiratory rate, tidal volume (TV) and positive end-expiratory pressure (PEEP) 

9 during the two types of scans. The data sets from the 4DCT scan were reconstructed in 

10 volumes at eight phases of the respiratory cycle. For this study we focus on the data 

11 sets from two of the phases, a volume near end expiration (EE) and a volume near end 

12 inspiration (EI). For the Xe-CT scan, 45 distinctive partial lung volumetric scans were 

13 performed at volume near end expiration, or the initial end expiration scan (EEq) to the 

14 last expiration scan (EE44). 

15 The nonlinear image registration is used to define the transformation Tl between the 

16 EE and EI in order to measure the regional lung ventilation from the 4DCT scan. The Xe- 

17 CT-based estimated regional lung ventilation is computed on the EEq by using Pulmonary 
IB Analysis Software Suite 11.0 (PASS) software by finding the constant of the exponential 

19 rise of the density from xenon gas wash- in over multiple breaths [TS]. The same nonlinear 

20 image registration is also applied to define the transformation T2 which maps the EEq 

21 to the EE so that the Xe-CT based estimate of ventilation can be mapped into the same 

22 coordinate system as that of the registration-based estimate of ventilation. Additional 

23 details on the registration algorithm and other processing steps are given below. 

24 2.2 Image Data Sets 

25 Appropriate animal ethics approval was obtained for these protocols from the University 

26 of Iowa Animal Care and Use Committee and the study adhered to NIH guidelines for 

27 animal experimentation. Four adult male sheep A, B, C, and D (with weights 44.0, 37.8, 
2B 40.4, and 46.7 kg) were used for this study. The sheep were anesthetized using intravenous 

29 pentobarbital and pancuronium to ensure adequate sedation and to prevent spontaneous 

30 breathing. Animals were positive pressure ventilated during experiments using a custom 

31 built dual Harvard apparatus piston ventilator designed for computer control. The 4DCT 

32 images were acquired with the animals in the supine position using the dynamic imaging 

33 protocol with a pitch of 0.1, slice coUimation of 0.6 mm, rotation time of 0.5 sec, slice 
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Figure 1: Figure shows the two types of images, a image pair of full lung volumetric phases 
EE and EI from a 4DCT scan and a Xe-CT scan acquired at end of expiration over 45 
respiratory cycles (EEq to EE44), which are analyzed during the processing. Transfor- 
mation Tl registers end inspiration (EI) to end expiration (EE) data and can be used to 
assess local lung function via calculations of three ventilation measures: specific air volume 
change by specific volume change (SAJ), specific air volume change by corrected Jacobian 
(SACJ), and specific air volume change by intensity (SAI). The 45 distinctive partial lung 
volumetric Xe-CT scans EEq to EE44 are used to calculate Xe-CT-based measure of spe- 
cific ventilation (sV). Transformation T2 maps the sV data into the coordinate system of 
the EE image (end expiration phase of the 4DCT scan) to allow direct comparison with 
the 4DCT and registration-based measures of ventilation. Both EE and EEq are at vol- 
umes near end inspiration. (Shaded boxes indicate CT image data; white boxes indicated 
derived or calculated data; thick arrows indicate image registration transformations being 
calculated; thin solid lines indicate other operations.) 
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thickness of 0.75 mm, slice increment of 0.5 mm, 120 kV, 400 mAs, and kernel B30f. The 
airway pressure signal was simultaneously recorded with the X-ray projections and images 
were reconstructed retrospectively at 0, 25, 50, 75, and 100% of the inspiration duration 
and 75, 50 and 25% of the expiration duration. The 0% (EE) and 100% (EI) inspiration 
phases were used for later ventilation measurements. A slab of twelve contiguous axial 
slices were imaged over 45 breaths for Xe-CT scans. Images were acquired using respiratory 
gating by triggering the scan during end-expiration with 80 keV energy (for higher density 
resolution, approximately 2 HU per 1% Xe), 160 mAs tube current, a 360° rotation, a 0.33 
sec scan time, and 2.4 mm slice thickness. Respiratory gating is achieved using a custom 
built Lab VIEW program which controls the ventilators and triggers the CT scanner. The 
respiratory rate (RR) for four animals ranged from 15 to 18 breaths per minute with an 
inspiratory-expiratory ratio of 1:1 which was sufficient to maintain a normocapnic state. 
Both of the two types of images were acquired with a matrix of 512 by 512 and without 
moving the animal between scans. 

2.3 Image Registration 

A tissue volume and vesselness measure preserving nonrigid registration (TVP) algo- 
rithm fL9[ [20] is used to estimate the transformations EI to EE and EEq to EE. The TVP 
algorithm minimizes the sum of squared tissue volume difference (SSTVD) [211 1221 [23l [2l] 
and vesselness measure difference (SSVMD), utilizing the rich image intensity information 
and natural anatomic landmarks provided by the vessels. This method has been shown to 
be effective at registering across lung CT images with high accuracy p!9l 120] . 

Let /i and I2 represent two 3D image volumes to be registered. The vector x = 
(xi, X2, x^Y' defines the voxel coordinate within an image. The algorithm finds the optimal 
transformation h that maps the template image Ji to the target image I2 by minimizing 
the cost function 

Ctotal = pI [y2(x)-Vl(h(x))]'cix + x / [F2(x)-Fi(h(x))]'cix. (1) 

where fl is the union domain of the lung regions in images Ii and 12- V2 and Vi are the 
tissue volumes as defined in Equation O F2 and Fi are the vesselness measures as defined 
in Equation El The transformation h is a (3 x 1) vector- valued function that maps a point 
h(x) in the target image to its corresponding location in the template image. The first 
integral of the cost function defines the SSTVD cost and the second integral of the cost 
function defines the SSVMD cost. 
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The SSTVD cost assumes that the measured Hounsfield units (HU) in the lung CT 
images is a function of tissue and air content. Following the air-tissue mixture model by 
Hoffman et al. [25], from the CT value of a given voxel, the tissue volume can be estimated 
as 

^(X) = -{^) rrjj^l7j = Kx)/?(/(x)), (2) 
U tissue ^ air 

and the air volume can be estimated as 

^'w = .(x) 2 = Kx)«(/(x)), (3) 

U tissue -n y air 

where i/(x) denotes the volume of voxel x and J(x) is the intensity of a voxel at position x. 
HUair and HUussue refer to the intensity of air and tissue, respectively. In this work, we 
assume that air is -1000 HU and tissue is HU. a(J(x)) = and /3(/(x)) = 

\ \ J/ n U tissue — nU air ' \ \ // 

HUr^~^-HU introduced for notational simplicity. Notice that a(/(x)) + /3(/(x)) = 1. 

Given (|2]), we can then define the SSTVD cost: 

CSSTVD = / [1^2(X)- Vi(h(x))]'dx (4) 

Jn 

= [ h(x)/3(/2(x))-j.i(h(x))/3(/i(h(x)))]^(ix, (5) 

The notation /i(h(x)) is interpreted as the image /i(x) deformed by the transformation 
h(x) and is computed using trilinear interpolation. The deformed volume element z/i(h(x)) 
is calculated using the Jacobian J(x) times the volume element i^2(x), i.e., i/i(h(x)) = 
J(x)z/2(x). 

Figure [2] shows an example of a cubic shaped region under deformation h from template 
image to target image. The region volumes are ui and i>2. The volumes can be decomposed 
into the tissue volume fraction and air volume fraction based on the mean voxel intensity 
within the cube. The small white sub volumes inside the cubes represent the tissue volume 
Vi and V2. Air volumes are represented by V{ and V2 (in blue). As the ratio of air to 
tissue decreases, the CT intensity of a voxel increases. The mean cube voxel intensities for 
the template, Ji, and target images, I2, are functions of the ratios of air to tissue volumes 
within the cubes. 

As the blood vessels branch to smaller and smaller diameters, the raw grayscale in- 
formation from vessel voxels provide very little contribution to guide the intensity-based 
registration. To better utilize the information of blood vessel locations, we use the vessel- 
ness measure based on the eigenvalues of the Hessian matrix of image intensity. Frangi's 
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Figure 2: Example of a region under deformation h(x) from template image to target 
image. Vi and V2 are tissue volumes in the regions. V{ and V^' are air volumes in the 
regions. Region volumes ui = Vi + V( and 1^2 = V2 + 



vesselness function [26j is defined as 



F(A) 
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if A2 < and A3 < 
otherwise 
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I A. 
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where Ra distinguishes between plate-like and tubular structures, Rb accounts for the 
deviation from a blob-like structure, and S differentiates between tubular structure and 
noise. The vesselness function has been previously widely used in vessel segmentations in 
lungs [271 [28] ^-iid in retinas [29j. a, /3, 7 control the sensitivity of the vesselness measure. 
The vesselness measure is rescaled to [0, 1] and can be considered as a probability-like 
estimate of vesselness features. For this study, a = 0.5, f3 = 0.5, and 7 = 5 and the 
weighting constants in the total cost were set as p = 1 and x = 0.2. These parameters are 
similar to those used in our previous work [191 12Q] . 

The transformation h(x) is a cubic B-splines transform: 

h(x)=x + 5^(/>,/3(3)(x-x,), (8) 



S 



Xl + Xl + A|, 



(7) 
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where 0j describes the displacements of the control nodes and /3'-^''(x) is a three-dimensional 
tensor product of basis functions of cubic B-Spline. A spatial multiresolution procedure 
from coarse to fine is used in the registration in order to improve speed, accuracy and 
robustness. The total cost in Equation [1] is optimized using a limited-memory, quasi- 
Newton minimization method with bounds (L-BFGS-B) [30] algorithm. The B-Splines 
coefficients are constrained so that the transformation maintains the topology using the 
sufficient conditions that guarantee the local injectivity of functions parameterized by 
uniform cubic B-Splines proposed by Choi and Lee |3T] . 

2.4 Regional Ventilation Measures from Image Registration 

After we obtain the optimal warping function h(x), we can calculate the regional ven- 
tilation, which is equal to the difference in local air volume change per unit time. The 
commonly-used ventilation measure is the specific ventilation sV which takes the initial 
air volume into account. The sV is equal to the specific air volume change sVol per unit 
time. Or in other words, in a unit time. 



Three different approaches for estimating are described below: 

Specific air volume change by specific volume change (SAJ): The SAJ regional 
ventilation measure is based on the assumption that local volume change is due to air 
volume change alone, and thus, any regional volume change is due only to local air volume 
change. Or, equivalently, the SAJ measure assumes that there is no tissue volume within 
the template or target volumes. Figure [3] illustrates such an assumption. Compared 
with the general condition in Figure El the region volume now is pure air volume, or 
equivalently, i^i = V{ and 1^2 = V2. In this case, the specific air volume change is equal 
to specific volume change. Since the Jacobian tells us the local volume expansion (or 
contraction), the regional ventilation can be measured by: 



Previously, SAJ has been used as an index of the regional function and was compared 
with Xe-CT estimates of regional lung function [1] . Regional lung expansion, as estimated 
from the Jacobian of the image registration transformations, was well correlated with 
xenon CT specific ventilation [H |6] (linear regression, average = 0.73). 



sV = sVol 




(9) 



SAJ 




J(x) - 1. 



(10) 
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Figure 3: Example of a given region under deformation h(x) from template image to 
target image, with the assumption of no tissue volume (Vi = V2 = 0). V{ and V2 are air 
volumes. 



1 Specific air volume change by corrected Jacobian (SAC J): Starting with (JTO] 

2 and expressing the air volumes V\'(h(x)) and V.^i'x) using the air-tissue mixture model ([2] 

3 and ([3]), we obtain the corrected Jacobian measure of region air volume change, SAC J, 



SACJ 



F/(h(x)) - ^^'(x) 



^/(h(x)) 
^'(x) 

z.i(h(x))a(/i(h(x))) ^ 
//2(x)a(/2(x)) 

4 As z/i(h(x)) = J(x)z/2(x), the specific air volume change is then 



(11) 
(12) 

(13) 



SACJ = J(x)^M^-l (14) 
a(/2(x)) 

= Ji^) ""';r-'f,'^f -i (15) 

nU tissue — 

5 If we assume that pure air is -1000 HU and pure tissue is HU, then specific air volume 
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1 change is 

SAC J = J(x)MM^-l. (16) 

^2(x) 

2 Compared to Equation [ini the term ^^j^l^^^^ is a correction factor that depends on the 

3 voxel intensities in the template and target images. The SACJ measure is illustrated in 

4 Figure [21 and represents the most general case of where there is both tissue volume and 

5 air volume change within the region. 



6 Specific air volume change by intensity change (SAI): The intensity-based mea- 

7 sure of regional air volume change SAI can be derived from the SACJ by assuming that 

8 tissue volume is preserved during deformation, or equivalently, that the tissue volume dif- 

9 ference A'\/(x) = Vi(h(x)) — V2(x) = 0. Under this assumption, Vi(h(x)) = V2(x) and we 

10 have 

i/i(h(x))/?(/i(h(x))) = z/2(x)/3(/2(x)), (17) 

11 and 

/?(/i(h(x))) 

12 Since z/i(h(x)) = J(x)zy2(x), with above equation, we have 

, , B(Ux.)) , , 



h{h{^))-HUai 

13 Substituting the above equation into equation [T51 yields 



(20) 



gj^j ^ /2(x) - HUair HUussue - /l(h(x)) _ ^ ^^^^ 
/l(h(x)) - HUair HUtissue " /2(x) 



(22) 



^ hi^) HUussue + HUg^hiH^)) - h{h{^)) HUussue ' HUg^hj^) 
(/l(h(x)) - HUa^r){HUussue " /2(x)) 

14 Finally, if we assume that pure air is -1000 HU and pure tissue is HU, then 

A(h(x)) -/2(x) , , 

SAI = 1000 , , ^ ' , 23 

/2(x)(/i(h(x)) + 1000) ^ ' 

15 which is exactly the result from Simon [32], Guerrero et al. [11], and Fuld et al. [33]. 

16 Figure H] illustrates the assumption with no tissue volume change in SAI. In Figure H] 
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Figure 4: Example of a given voxel under deformation h(x) from template image to target 
image, with the assumption of no tissue volume change. Notice the tissue volume Vi = V2 
under this assumption. V( and V2 are air volumes. 

as the region volume changes from ui to z/2, the tissue volume inside the cube remains the 
same (Vi = V2). 

Difference of specific air volume change (DSA) and difference of tissue vol- 
ume (DT): To investigate the relationship between the measurements of specific air 
volume changes and the tissue volume change, we can calculate the difference between 
equation (fT5l) and equation (!22|) and define the difference of specific air volume change 
(DSA) between SACJ and SAI, and the difference of tissue volume (DT) as: 



DSA= \SACJ-SAI\ 



(24) 



DT = |ri(h(x))-y2(x)| (25) 

= |z.i(h(x))/?(/i(h(x)))-/.2(x)/?(/2(x))| (26) 

= |J(x)z.2(x)/?(/i(h(x)))-/y2(x)/?(/2(x))| (27) 

I J(x)(/i(h(x)) - HU,,r) - (/2(X) - HU,,r) , .... 

^ I 2lXj — —-zr- I (28 j 

^ ^tissue ^ ^ air 

Again, if we assume that air is -1000 HU and tissue is HU, then the tissue volume 
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(29) 



2.5 Computational Setup 

Processing starts by identifying the lung regions in all images using the Pulmonary Work- 
station 2.0 (VIDA Diagnostics, Inc., Iowa City, lA). The Xe-CT estimate of sV is computed 
in the coordinates of the EEq using Pulmonary Analysis Software Suite 11.0 (PASS) |T8] at 
the original image size of 0.5 mm x 0.5 mm x 2.4 mm voxels. Overlapping 1x8 regions 
of interest (ROI) are defined in the lung region on each 2D slice. All images, including the 
EE, EI, EEq and their corresponding lung region masks or sV map, are then resampled to 
a voxel size of 1 mm x 1 mm x 1 mm. After preprocessing, EI is registered to EE using 
the TVP registration for measuring the regional ventilation from these two phases in a 
4DCT scan. The resulting transformation is used to estimate the SAJ, SACJ and SAL 
Then EEq is registered to EE using TVP registration as well to map the sV to the same 
coordinate system as that of the SAJ, SACJ and SAL Due to the fact that the denominator 
of equation ( IT5|) and (!22|) may become zero, we eliminated from consideration any points 
that have the absolute value of the denominator become less than 0.001. For the TVP 
registration, the multiresolution strategy is used in the processing and it proceeds from 
low to high image resolution starting at one-eighth the spatial resolution and increases by 
a factor of two until the full resolution is reached. Meanwhile, a hierarchy of B-spline grid 
spaces from large to small is used. The finest B-spline grid space used in the experiments 
is 4 mm. The images and image grid space are refined alternatively. 

2.6 Assessment of image registration accuracy 

A semi-automatic landmark system is used for landmark selection and matching [34] . This 
system first uses an automatic landmark detection algorithm to find the landmarks in the 
EE image. The algorithm automatically detects "distinctive" points in the target image 
as the landmarks based on a distinctiveness value D{p). Around each point p, 45 points, 
gi, . . . ,^45 are uniformly distributed on a spherical surface. A region of interest ROI{qi) 
is compared with the corresponding region of interest ROI{p) around the original point, 
and then combined with its gradient value to calculate the distinctiveness value D{p). 

The same system is then applied to guide the observer to match landmarks in the target 
image with corresponding landmarks in the template image. Each landmark-pair manually 
annotated by the observer is added to a thin-plate-spline to warp the template image. The 
system utilizes the warped image to estimate where the anatomic match will be located 
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for a new landmark point presented to the observer, therefore the observer can start the 
matching from a system estimated location. Thus, as the warped image becomes more 
accurate by the new added landmarks, the task of the observer becomes easier. 

For each animal, after 200 anatomic landmarks are identified in the EE image, the 
observer marks the locations of the voxels corresponding to the anatomic locations of the 
landmarks in the EI image. For each landmark, the actual landmark position is compared 
to the registration-derived estimate of landmark position and the error is calculated. With 
the evaluated accuracy of transformation from the lung image registration algorithm, the 
resulting regional ventilation measures estimated using the transformation can be then 
compared to Xe-CT estimated ventilation. 

2.7 Compare Registration Regional Ventilation 
Measures to Xe-CT Estimated Ventilation 

In our previous work |H |6], regional lung expansion, as estimated from the Jacobian of 
the image registration transformations, was compared with Xe-CT-based sV. The analysis 
was conducted by evaluating Jacobian value between a pair of lung volumes from static 
(multiple airway pressures) and dynamic image data sets, and comparing the Jacobian 
along the y (ventral-dorsal) axis. While the correlation between the Jacobian value and 
sV reflect the fact that regional expansion estimated from image registration can be used 
as an index as regional lung function, the spatial resolution of the analysis method em- 
ployed were likely not sufficient to distinguish the differences between regional ventilation 
measures as we have described in Section 12.41 Therefore, to locally compare the regional 
ventilation measures, the corresponding region of Xe-CT image EEq in the EE is divided 
into approximately 100 non-overlapping cubes with size of 20 mm x 20 mm x 20 mm. 
We compare the average regional ventilation measures (SAJ, SACJ and SAI) to the corre- 
sponding average sV measurement from Xe-CT images within each cube. The correlation 
coefficients between any two estimates (SAJ-sV, SACJ-sV or SAI-sV) are calculated by 
linear regression. To compare two correlation coefficients, the Fisher Z-transform of the r 
values is used and the level of significance is determined [35]. The relationship between the 
specific air volume change and difference of tissue volume is also studied in four animals 
by linear regression analysis. 
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3 Results 



3.1 Registration Accuracy 

For each animal, approximately 200 automatically identified landmarks within the lungs 
are used to compute registration accuracy. The landmarks are widely distributed through- 
out the lung regions. Figure shows an example of the distribution of the landmarks in 
animal D for both the EE and EI images. The coordinate of each landmark location is 
recorded for each image data set before and after registration for all four animals. Fig- 
ure |6] shows the landmark distance before and after registration for four animals. The grey 
boxes show the magnitude of respiratory motion during the tidal breathing. For all four 
animals, before registration, the average landmark distance is 6.6 mm with a minimum 
distance of 1.0 mm, maximum distance of 14.6 mm, and standard deviation of 2.42 mm. 
After registration, the average landmark distance is 0.4 mm with a minimum distance of 
0.1 mm, a maximum distance of 1.6 mm, and a standard deviation 0.29 mm. The trends 
for all animals are consistent and the results demonstrate that the registrations produced 
good anatomic correspondences. All registrations were examined and it was confirmed 
that all Jacobian values had positive values. 

Figure |5t^a) shows the location of the EEq (Xe-CT) slab overlaid on the EE image. 
Figure [7] shows an example of the image registration result from the EEq image to the EE 
image. The first row shows the misalignment between the images before image registration. 
Though the images were acquired without moving the animal between the scans, there is 
still non-rigid deformation between scans as shown in Fig[71^d), as the black and white 
regions represent the large intensity difference between Fig. [Tl^a) and (b). In addition, the 
slice thicknesses were quite different which causes partial volume artifacts. After image 
registration, the EEq image is aligned to the EE image, and the resulting difference image 
(shown in Fig.[7](e)) is near zero. The transformation from the EEq to the EE image allows 
us to map the Xe-CT sV into the coordinate system of EE image. Note that since the 
regions outside the lung are not included in the registration process, the mediastinum and 
other body tissues are not aligned. Also note that the dorsal region of the lung shows 
a intensity difference after registration. This is due mainly to the gradual progression of 
atelectasis and tissue edema during the course of the experiment. 
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(a) 



(b) 



Figure 5: 3D view of the landmarks in: (a) EE with EEq and (b) EI. The dark region 
below the carina in (a) is the EEq and the spheres are the automatically defined landmarks. 

1 3.2 Registration Estimated Ventilation Compared 

2 to Xe-CT Estimated Ventilation 

3 Figure IH!^a) shows a comparison between the registration-derived indices of ventilation 

4 and the Xe-CT estimated sV in cube-shaped regions of interest for animal D. The corre- 

5 sponding Xe-CT regions in the EE are divided into about 100 cubes. Figure Mjo) is the 

6 Xe-CT estimate of sV. Figure [8]^c), (d), (e) are the corresponding registration ventilation 

7 measures SAJ, SACJ and SAL The regions with edema are excluded from the compari- 

8 son. Figure El^b) to (d) all show noticeable similar gradient in the ventral-dorsal direction. 

9 Notice that the color scales are different in each map and are set based on the range of 

10 values from the appropriate plot in Fig. |9l 

11 Figure [9] shows scatter plots comparing the registration ventilation measures and the Xe- 

12 CT ventilation sV in all four animals. The SACJ column shows the strongest correlation 

13 with the sV (average = 0.82). The SAJ, which is directly related to Jacobian as 

14 SAJ = J — 1, also shows good correlation with the sV (average = 0.75). The intensity- 
is based measure SAI shows the lowest correlation with the sV (average = 0.58). 

16 Table [1] shows the results of comparing the r values from SACJ vs. sV and SAI vs. sV. 

17 All four animals show that the correlation coefficient from SACJ vs. sV is significantly 
IB stronger than it from SAI vs. sV. Similarly, table |2] shows the results of comparing the 
19 r values from SAJ vs. sV and SAI vs. sV. The registration ventilation measure SAJ also 
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□ Before T1 Registration 

□ After T1 Registration 
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Figure 6: Landmarks distances of the registration pair EI to EE for all four animals. 
Boxplot lower extreme is first quartile, boxplot upper extreme is third quartile. Median 
is shown with solid horizontal line. Whiskers show either the minimum (maximum) value 
or extend 1.5 times the first to third quartile range beyond the lower (upper) extreme of 
the box, whichever is smaller (larger). Outliers are marked with circles. 



19 



(a) 




(d) (e) 

Figure 7: Visualization of the result of the transformation that maps the Xe-CT estimated 
ventilation sV to the EE coordinate system: (a) EE, (b) EEq, (c) deformed EEo after 
registration, (d) intensity difference between EE and EEq before registration, (e) intensity 
difference between EE and EEq after registration. 
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shows a significantly stronger correlation with sV than SAL Comparing the r values from 
SACJ vs. sV and SAJ vs. sV, only animal B and C show that the SACJ has significantly 
stronger correlation with sV than SAJ. 

To analyze the effect of the size of the region of interest, the corresponding region 
of Xe-CT image EEq in the EE is divided into about 30 slabs along the ventral-dorsal 
direction with size of 150 mm x 8 mm x 40 mm as similarly in our previous work |1[ [6]. 
Figure [10] shows the scatter plots between the registration ventilation measures and the 
Xe-CT ventilation sV similar as Fig. [10] but in larger ROIs. The SACJ column shows 
the strongest correlation with the sV (average = 0.92). Both the SAJ and SAI show 
good correlation with sV as well (average = 0.88 and = 0.87). However, though 
the average value still show the SACJ has the highest correlation with Xe-CT based 
sV, Table [1] and [5] show that with larger averaging region as defined slabs, there is no 
significant difference between the correlation coefficients from SACJ vs. sV and SAI vs. 
sV, between SAJ vs. sV and SAI vs. sV, or SACJ vs. sV and SAJ vs. sV as in Table [6] 

Figure [11] shows the scatter plots between DSA (the absolute difference of the value 
between the SACJ and SAI) and the DT (the absolute difference of the tissue volume) 
with linear regression in all four animals (average = 0.86). From the equation f lTSj) 
and ( l22l) . we know that the SAI measurement assumes no tissue volume change in a given 
region being registered, which may not be valid assumption in all lung regions. Figure [TT] 
shows that as the tissue volume change increases, the difference between the measures of 
regional ventilation from SACJ and SAI increases linearly in all four animals. It indicates 
that the lower correlation of SAI with sV compared with SACJ with sV may be caused 
by the tissue volume change between two volumes. 

4 Summary and Conclusions 

We described three measures to estimate regional lung tissue ventilation from tissue volume 
and vesselness preserving image registration of CT images. These measures have been 
compared with each other, and compared to Xe-CT estimates of specific ventilation. We 
examined the assumption of constant tissue volume between registered lung regions, and 
demonstrated that the difference between two of the registration derived measures (SACJ 
and SAI) may be explained by differences in tissue volume between the lung regions being 
compared with registration. 

The tissue volume and vesselness preserving non-linear registration algorithm was used 
to match the EI image to the EE image to produce the registration deformation field and 
estimates of regional ventilation. It was used to register the EEq image to the EE image for 
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(a) 




0.0 1.0 2.0 3.0 4.0 5.0 6.0 0.0 1.2 2.3 3.5 4.7 5.9 7.0 



(b) (c) 
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(d) (e) 

Figure 8: Comparison of the regional ventilation measures for animal D. (a): EE with 
color coded cubes showing the sample region, (b), (c), (d) and (e): color map of the sV, 
SAJ, SAC J and SAL Note that the color scales are different for (b)-(e), and are set based 
on the range of values from the appropriate plot in Fig. [9l The results were similar for 
the other three animals 22 
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Figure 9: Small cube ROls with size 20 mm x 20 mm x 20 mm results for registration 
estimated ventilation measures compared to the Xe-CT estimated ventilation sV in scatter 
plot with linear regression in four animals. The first column is the SAJ vs. sV. The second 
column is the SACJ vs. sV. The third column is the SAI vs. sV. 
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Figure 10: Large slab ROls with size 150 mm x 8 mm x 40 mm results for registration 
estimated ventilation measures compared to the Xe-CT estimated ventilation sV in scatter 
plot with hnear regression in four animals. The first column is the SAJ vs. sV. The second 
column is the SACJ vs. sV. The third column is the SAI vs. sV. 
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Animal 


Correlation 


Correlation with sV 


Number of 


SACJ vs. SAI 




pair 


(r value) 


samples 


•p value 


A 


SACJ vs. sV 


0.88 


83 


p 0.0001 




OAT ^ ^\ T 

bAl vs. sV 


U.DO 






B 


SACJ vs. sV 


0.93 


119 


p <= 1.18e-6 




SAI vs. sV 


0.77 






C 


SACJ vs. sV 


0.89 


86 


p <= 0.0075 




SAI vs. sV 


0.78 






D 


SACJ vs. sV 


0.92 


110 


p <= 0.0017 




SAI vs. sV 


0.83 







Table 1: Comparison of ventilation measures between SACJ and SAI in small cube ROIs 
with size 20 mm x 20 mm x 20 mm. 



Animal 


Correlation 


Correlation with sV 


Number of 


SAJ vs. SAI 




pair 


(r value) 


samples 


p value 


A 


SAJ vs. sV 


0.86 


83 


p <= 0.0005 




SAI vs. sV 


0.65 






B 


SAJ vs. sV 


0.89 


119 


p <= 0.002 




SAI vs. sV 


0.77 






C 


SAJ vs. sV 


0.78 


86 


p <= 0.5 




SAI vs. sV 


0.78 






D 


SAJ vs. sV 


0.92 


110 


p <= 0.0017 




SAI vs. sV 


0.83 







Table 2: Comparison of ventilation measures between SAJ and SAI in small cube ROIs 
with size 20 mm x 20 mm x 20 mm. 
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Animal 


Correlation 


Correlation with sV 


Number of 


SACJ vs. SAJ 




pair 


(r value) 


samples 


p value 


A 


SACJ vs. sV 


0.88 


83 


p <= 0.302 




OAT ^ ^\ T 

bAJ vs. sV 


0.8d 






B 


SACJ vs. sV 


0.93 


119 


V <= 0.035 




SAJ vs. sV 


0.89 






C 


SACJ vs. sV 


0.89 


86 


p <= 0.007 




SAJ vs. sV 


0.78 






D 


SACJ vs. sV 


0.92 


110 


p <= 0.5 




SAJ vs. sV 


0.92 







Table 3: Comparison of ventilation measures between SACJ and SAJ in small cube ROIs 
with size 20 mm x 20 mm x 20 mm. 



Animal 


Correlation 
pair 


Correlation with sV 
(r value) 


Number of 
samples 


SACJ vs. SAI 
p value 


A 


SACJ vs. sV 


0.95 


17 


p <= 0.5 




SAI vs. sV 


0.95 






B 


SACJ vs. sV 


0.99 


23 


p <= 0.005 




SAI vs. sV 


0.95 






C 


SACJ vs. sV 


0.94 


23 


p <= 0.15 




SAI vs. sV 


0.89 






D 


SACJ vs. sV 


0.95 


25 


p <= 0.28 




SAI vs. sV 


0.93 







Table 4: Comparison of ventilation measures between SACJ and SAI in large slab ROIs 
with size 150 mm x 8 mm x 40 mm. 
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Animal 


Correlation 


Correlation with sV 


Number of 


SAJ vs. SAI 




pair 


(r value) 


samples 


p value 


A 


SAJ vs. sV 


0.95 


17 


p <= 0.5 




OAT ^ ^\ T 

bAl vs. sV 


0.95 






B 


SAJ vs. sV 


0.99 


23 


p <= 0.005 




SAI vs. sV 


0.95 






C 


SAJ vs. sV 


0.94 


23 


j9 <= 0.16 




SAI vs. sV 


0.89 






D 


SAJ vs. sV 


0.95 


25 


p <= 0.28 




SAI vs. sV 


0.93 







Table 5: Comparison of ventilation measures between SAJ and SAI in large slab ROIs 
with size 150 mm x 8 mm x 40 mm. 



Animal 


Correlation 


Correlation with sV 


Number of 


SACJ vs. SAI 




pair 


(r value) 


samples 


•p value 


A 


SACJ vs. sV 


0.95 


17 


p <= 0.5 




SAJ vs. sV 


0.95 






B 


SACJ vs. sV 


0.99 


23 


p <= 0.5 




SAJ vs. sV 


0.99 






C 


SACJ vs. sV 


0.94 


23 


p <= 0.5 




SAJ vs. sV 


0.94 






D 


SACJ vs. sV 


0.95 


25 


p <= 0.5 




SAJ vs. sV 


0.95 







Table 6: Comparison of ventilation measures between SACJ and SAJ in large slab ROIs 
with size 150 mm x 8 mm x 40 mm. 



27 



0.7 

0.6 
§0.5 

-I—' 

5.0.4 

^0.3 
< 

WD 0.2 
Q 

0.1 
0.0 



y = 4.845x + 0.002 
= 0.88 




0.7 



0.00 002 OM 006 O08 1 



DT (1 /inspiration) 

(a) 



y = 5.056x-0.001 



= 0.893 




0.00 0.02 0.04 0.06 0.08 0.10 



DT (1/inspiration) 

(b) 



0.7 

0.6 

gO.5 

5.0.4 
w 

^0.3 
< 

CO 0.2 
Q 

0.1 
0.0 



y = 4.284x + 0.004 
= 0.889 




0.00 0.02 0.04 0.06 0.08 
DT (1/inspiration) 



0.10 



0.7 
0.6 



y = 4.387x + 0.008 
= 0.839 




0.00 0.02 0.04 0.06 0.08 
DT (1/inspiration) 

(d) 



0.10 



Figure 11: Linear regression analysis between DSA and DT. (a) to (d): DSA (the absolute 
difference of the value between the SACJ and SAI) compared to DT (the absolute difference 
of the tissue volume) in animals A, B, C and D. 
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comparing the three ventilation measures to the Xe-CT based sV. About 200 anatomical 
landmarks were identified and annotated in each data set to evaluate registration accuracy. 
The average landmark error is on the order of 1 mm after registration. 

The ventilation measures SAJ, SACJ, and SAI were derived using a simple model of a 
lung region, containing a mixture of air and tissue, which deforms during inspiration or 
exhalation. The SAJ measure, which is a linear function of the Jacobian of the registration 
displacement field, measures regional ventilation based on the assumption that the lung 
region contains only air (i.e., no tissue volume). The SACJ is the most general form of 
the three measures and is based on model where both the air and tissue volumes can 
change during inspiration and exhalation. Finally, the SAI measure is computed based on 
intensity change alone, and assumes that the region may have a tissue volume that is non- 
zero, but this volume does not change during inspiration and exhalation. Thus, the SAJ 
measure relies solely on the volume change information computed from the Jacobian of 
the deformation field and the SAI measure relies solely on the change in region intensity as 
measured by the CT. The SACJ measure uses a more complete model of local ventilation, 
and combines the geometric information from the Jacobian with the density information 
calculated from the change in region intensity. 

The three registration-based ventilation measures and the Xe-CT sV measurement were 
averaged and compared in cubic-shaped regions of interest. In 20 mm x 20 mm x 20 
mm ROIs, the SACJ shows significantly higher correlation with Xe-CT sV than the SAI 
in all four animals. By studying the difference between the SACJ and SAI measures and 
the tissue volume difference estimated by the CT intensity change, we showed that the 
difference between SACJ and SAI may be explained by the constant tissue volume as- 
sumption implicit in the SAI model f|T9|) . From Fig. [Ill we see that the difference between 
the SACJ and SAI measures is approximately linearly-related to the estimated tissue vol- 
ume change. Tables [1] and [2] show that the both the SACJ and SAJ have significantly 
better correlation with sV than the SAI. This is consistent with the findings by Kabus 
et al. [36] who showed that the Jacobian-based measure of ventilation has less error than 
the intensity-based ventilation measure, using the segmented total lung volume as a global 
comparison. Though all the regional ventilation measures and Xe-CT based sV from the 
linear regression analysis in Fig. [TH] show high correlations. Tables [H and [5] show that 
there is no significant difference in the correlation with sV between the Jacobian-based 
measures and intensity-based measure. This result indicates that the validation methods 
using global comparison such as segmented total volume may not be able to distinguish 
the Jacobian-based measure and the intensity-based measure. 

The comparison of the ventilation measures was limited to the resolution of 20 mm x 
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20 mm X ROIs. As the size of the ROIs decreases, the correlation between the ventilation 
measures with Xe-CT based sV decreases. This may be due to the underlying noise of 
the Xe-CT measurement of ventilation or the decreased sensitivity of registration based 
measure to local ventilation heterogeneity which is relative to the case. Additional Xe-CT 
image analysis work including using multi-compartment models, thinner slice, and inter- 
phase registration to improve sV measurement are required to reduce the noise in Xe-CT 
based sV measurement. 

To compare with the intensity-based ventilation measure used in previous work in Si- 
mon [52], Guerrero et al. [TT], and Fuld et al. [33], we followed the assumption that HUair 
is -1000 HU and HUussue is HU (equaling water [25]) in this work. The ventilation mea- 
sures were calculated under the assumption that HUair is -1000 HU and HUussue is 55 HU, 
which are the values used by Yin et al. in [22]. Our analysis shows that the correlation 
coefficients between any two estimates (SAJ-sV, SACJ-sV or SAI-sV) change less than 
1% with two different HUussue values. However, it would be important to have sensitivity 
analysis in the future to compare different ventilation measures against intensity changes. 

The image registration algorithm used to find the transformation from EI to EE for 
measurement of regional ventilation produces accurate registrations by minimizing the 
tissue volume and vesselness measure difference between the template image and the target 
image. It would be interesting to compare different image registration algorithms and their 
effects on the registration-based ventilation measures. For example, if two registration 
algorithms achieve the similar landmark accuracy, the one does not preserve tissue volume 
change may show even larger difference in the SACJ and SAI measures than the results 
using TVP as described above. 

In conclusion, with the same deformation field by the same image registration algorithm, 
a significant difference between the Jacobian/-corrected Jacobian-based ventilation mea- 
sures and the intensity-based ventilation measure is found in a regional level using Xe-CT 
based ventilation measure sV. The ventilation measure by corrected Jacobian SACJ gives 
best correlation with Xe-CT based sV and the correlation is significantly higher than from 
the ventilation by intensity SAI indicating the ventilation measure by corrected Jacobian 
SACJ may be a better measure of regional lung ventilation from image registration of 
4DCT images. 
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